R Notebook

https://jakubnowosad.com/posts/2024-10-20-spatcomp-bp2/

library(terra)

terra 1.7.23

ndvi2018_tartu = rast("https://github.com/Nowosad/comparing-spatial-patterns-2024/raw/refs/heads/main/data/ndvi2018_tartu.tif")
ndvi2023_tartu = rast("https://github.com/Nowosad/comparing-spatial-patterns-2024/raw/refs/heads/main/data/ndvi2023_tartu.tif")
plot(ndvi2018_tartu, main = "Tartu (2000)")

plot(ndvi2023_tartu, main = "Tartu (2018)")

ndvi_diff = ndvi2023_tartu - ndvi2018_tartu
terra::plot(ndvi_diff)

plot_div = function(r, ...){
  r_range = range(values(r), na.rm = TRUE, finite = TRUE)
  max_abs = max(abs(r_range))
  new_range = c(-max_abs, max_abs)
  plot(r, col = hcl.colors(100, palette = "prgn"), range = new_range, ...)
}
plot_div(ndvi_diff)

ndvi_diff = ndvi2023_tartu - ndvi2018_tartu
hist(ndvi_diff)

library(yardstick)
ndvi_rmse = rmse_vec(values(ndvi2023_tartu)[, 1], values(ndvi2018_tartu)[, 1])
ndvi_rmse

[1] 0.2191853

library(diffeR)

Loading required package: ggplot2

ndvi_mad = MAD(ndvi2023_tartu, ndvi2018_tartu)

[1] "The mean of grid1 is less than the mean of grid2"

ndvi_mad$Total

[1] 0.184306

ndvi_cor = focalPairs(c(ndvi2023_tartu, ndvi2018_tartu), w = 5,
                      fun = "pearson", na.rm = TRUE)
plot_div(ndvi_cor)

library(geodiv)

Attaching package: ‘geodiv’

The following objects are masked from ‘package:terra’:

    rotate, sds

window = matrix(1, nrow = 5, ncol = 5)
ndvi2018_tartu_sa_mw = focal_metrics(ndvi2018_tartu, window = window,
                               metric = "sa", progress = FALSE)
ndvi2023_tartu_sa_mw = focal_metrics(ndvi2023_tartu, window = window,
                               metric = "sa", progress = FALSE)
ndvi_diff_sa_mw = ndvi2023_tartu_sa_mw$sa - ndvi2018_tartu_sa_mw$sa
plot_div(ndvi_diff_sa_mw)

plot(ndvi2023_tartu)

plot(ndvi2023_tartu_sa_mw$sa)

library(rasterdiv)
ndvi2018_tartu_int = as.int(ndvi2018_tartu * 100)
ndvi2023_tartu_int = as.int(ndvi2023_tartu * 100)
ndvi2018_tartu_rao = paRao(ndvi2018_tartu_int, window = 5, progBar = FALSE)

Warning: Simplify=0. Rounding data to 0 decimal places.

Processing alpha: 1 Moving Window: 5


Processing alpha: 1 Moving Window: 5
Error in utils::combn(p, m = 2, FUN = prod, na.rm = TRUE) : n < m

ndvi_diff = ndvi2023_tartu - ndvi2018_tartu
ndvi_diff_autocor = autocor(ndvi_diff, method = "moran", global = FALSE)
plot_div(ndvi_diff, main = "Difference")

plot_div(ndvi_diff_autocor, main = "Moran's I of the difference")

library(SSIMmap)
ndvi_ssim = ssim_raster(ndvi2018_tartu, ndvi2023_tartu, global = FALSE, w = 3)
plot_div(ndvi_ssim[[1]])

library(waywiser)
cell_sizes = seq(50, 300, by = 50)
ndvi_multi_scale = ww_multi_scale(truth = ndvi2018_tartu, estimate = ndvi2023_tartu,
                                 metrics = list(yardstick::rmse), 
                                 cellsize = cell_sizes,
                                 progress = FALSE)

ℹ The package "exactextractr" is required.
✖ Would you like to install it?

1: Yes
2: No

1

Installing package into ‘C:/Users/stahlm/AppData/Local/R/win-library/4.2’
(as ‘lib’ is unspecified)
trying URL 'https://cran.rstudio.com/bin/windows/contrib/4.2/exactextractr_0.10.0.zip'
Content type 'application/zip' length 1448767 bytes (1.4 MB)
downloaded 1.4 MB

package ‘exactextractr’ successfully unpacked and MD5 sums checked

The downloaded binary packages are in
    C:\Users\stahlm\AppData\Local\Temp\RtmpWCibE9\downloaded_packages

Warning: No CRS specified for polygons; assuming they have the same CRS as the raster.Warning: No CRS specified for polygons; assuming they have the same CRS as the raster.Warning: No CRS specified for polygons; assuming they have the same CRS as the raster.Warning: No CRS specified for polygons; assuming they have the same CRS as the raster.Warning: No CRS specified for polygons; assuming they have the same CRS as the raster.Warning: No CRS specified for polygons; assuming they have the same CRS as the raster.Warning: No CRS specified for polygons; assuming they have the same CRS as the raster.Warning: No CRS specified for polygons; assuming they have the same CRS as the raster.Warning: No CRS specified for polygons; assuming they have the same CRS as the raster.Warning: No CRS specified for polygons; assuming they have the same CRS as the raster.Warning: No CRS specified for polygons; assuming they have the same CRS as the raster.Warning: No CRS specified for polygons; assuming they have the same CRS as the raster.

ndvi_multi_scale

library(diffeR)
ndvi_mad = MAD(ndvi2023_tartu, ndvi2018_tartu, eval = "multiple")

[1] "The mean of grid1 is less than the mean of grid2"

ndvi_mad

library(SSIMmap)
ssim_raster(ndvi2018_tartu, ndvi2023_tartu, global = TRUE)

SSIM: 0.63845 SIM: 0.90432 SIV: 0.90778 SIP: 0.75671